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ABSTRACT 

We investigate the use of a logarithmic density variable in estimating the Lagrangian displacement 
field, motivated by the success of a logarithmic transformation in restoring information to the matter 
power spectrum. The logarithmic relation is an extension of the linear relation, motivated by the 
continuity equation, in which the density field is assumed to be proportional to the divergence of 
the displacement field; we compare the linear and logarithmic relations by measuring both of these 
fields directly in a cosmological A^-body simulation. The relative success of the logarithmic and 
linear relations depends on the scale at which the density field is smoothed. Thus we explore several 
ways of measuring the density field, including Cloud-In-Cell smoothing, adaptive smoothing, and 
the (scale-independent) Delaunay tessellation, and we use both a Fourier space and a geometrical 
tessellation approach to measuring the divergence. We find that the relation between the divergence 
of the displacement field and the density is significantly tighter with a logarithmic density variable, 
especially at low redshifts and for very small (^2 h^^Mpc) smoothing scales. We find that the grid- 
based methods are more reliable than the tessellation-based method of calculating both the density 
and the divergence fields, though in both cases the logarithmic relation works better in the appropriate 
regime, which corresponds to nonlinear scales for the grid-based methods and low densities for the 
tessellation-based method. 

Subject headings: cosmology: theory — large scale structure of Universe — methods; statistical 



1. INTRODUCTION 

Under the influence of gravity, the initially-small fluc- 
tuations in the primordial density field grow to become 
the hierarchical and nonlinear structures observed to- 
day. We can study the late-time clustering of matter 
by evolving a density field described by an initial power 
spectrum and comparing the simulated structures to ob- 
served galaxies in a statistical way. However, if we were 
able to undo the smearing effects of nonlinear evolu- 
tion and reverse the formation of structure, we would be 
able to reconstruct the primordial power spectrum from 
an observed distribution of galaxi es. Many reconstruc- 
tion methods have been developed (Hamilton et al .|1991 



Weinberg 1992; Peacock & Dodds 1994, 1996; I-risch et 
al. 2002; B renier et al. 2003; Mohayaec et al. 2006),Irom 
simple linear reconstruction to mor e complex methods, 
with varying degrees of success (see ( Narayanan & Croft 
19991 for a comparative review). For exampl e, Monge- 
A mpere-Kantorovitch (MAK) reconstruction ( Frisc h et 



ari|2002l iBrenier et all 20031 iMohayaee et al.pooei was 



successfully able to r eproduce the observ ed velocity field 
of the Local Group (Lavaux et al. 2010). A reconstruc- 
tion method that uses the linear Zei'dovich approxima- 
tion of the Lagrangian displacement has also been pro- 
posed to enhance t he measurement of baryon acoustic 
oscillations (BAOs; Eisenstein et al. 2007| |No h et "aL 
20091 iPadmanabhan et al.||i^00j)[ |WMte| |!^blO; Me hta et, 
al.||201ip . The displacement field is a crucial quantity in 
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most of these methods, and is of interest in its own right. 

There is reason to believe that a logarithmic trans- 
formation of the density field may aid estimation of the 
displacement field. The mass distribution has been suc- 
cessfully described by a lognorm al field that evolve d from 
Gaussian initial conditions (Coles fc Jones 1991), with 
some evidenc e for a skewed lognormal field in the non- 
linear regime ( Colombi|1994 ). More recently, it has been 
found that the power spectrum of the log-transformed 
density field contains more Fisher information than the 
usual p ower spectrum at sma ll scales by up to a factor 



of ^10 (Neyrinck et al. 2009). A modified logarithmic 



transform has also been shown to increase the precision 
of the power sp ectrum of the non linear weak lensing con- 
vergence field (Seo et al. 20101. The log-transformed 



density field is also more effective in constraining cosmo 
logical parameters than the standar d dens ity field when 
using the power spectrum (Neyrinck||2011 ). 

This paper investigates the effect of a logarithmic 
transform of the density field on the relation between 
the density and displacement fields. We do this by mea- 
suring the divergence of the displacement field, which in 
linear theory is proportional to the negative density con- 
trast 5, using both Eulerian and Lagrangian techniques 
in a cosmological TV-body simulation. We compare the 
linear and logarithmic approaches and evaluate their de- 
pendence on redshift and smoothing scale, both of which 
affect the applicability of linear theory. The linear and 
logarithmic approximations for the displacement field are 
derived in Section [2] and three methods of measuring the 
divergence of the displacement field and the density field 
are presented in Section [3] Concluding remarks are given 
in Section m 
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2. THEORY 

In this section, we compare the hnear approximation 
of the continuity equation to a logarithmic approxima- 
tion. The large-scale dynamical evolution o f structure 
is des cribed by the Zel'dovich approximation ( Zel'dovich 
1970| , 

* = q ZJV^, 



X = q 



(1) 



where x is the comoving Eulerian coordinate, q the La- 
grangian coordinate, ^ the displacement, D the lin- 
ear growth function, and $ the gravitational poten- 
tial. The linear Poisson equation relates the density con- 
trast, S = (p — p) / p, to the gravitational potential by 
6 = £) V^$, thus the divergence of the linear Zel'dovich 
displacement can be related to the density contrast: 



V * 



(2) 



The time derivative of the Zel'dovich approximation in 
expanding coordinates (r = ax, where r is the physi- 
cal coordinate) allows us to relate this expression to the 
continuity equation. We first note that the derivative 
consists of two components, r = dx -I- ax, where the first 
component describes the rate of expansion and the sec- 
ond the peculiar velocity of galaxies (or mass particles): 
V = ax = a' p. The conti nuity equation can be written 
as (see, e.g., |Peebles]|1980[ p 48) 



djl + S) 
dt 



+ -V- [(l + <5)v] 



0. 



(3) 



-V • V = 0, 



(4) 



To first in 6 order this becomes 

m 1 
Ik 

thus we have 5 ~ —V - 4/' and have recovered EquationjSj 

The continuity equation (Eqn. |3| may also be simpli- 
fied by keeping the (1 -f 5) together, resulting in a loga- 
rithmic derivative of (1 + 5): 



1 d{l + S) 
(1 + 5) dt 



-V • V = 0. 



(5) 



This gives a logarithmic expression for the divergence of 
the Lagrangian displacement. 



V • * = - In (1 + 5) + C, 



(6) 



where the C term appears because ln(l -f 5) is not a 
zero-mean-field, so we must take into account a non-zero 
integration constant. We find that this constant is well- 
approximated by the rather natural value of (ln(l -I- 5)). 
We now may compare the linear relation, given by Eqn. [21 
to the logarithmic relation, given by Eqn. [6j 

3. RESULTS 

We measure the density field and V • * in a 
200 h-^Mpc, 256^ particle Gadget ( ISpringel et al||2001 1 
cold dark matter simulation with standard ACUM cos- 
mology {VLm = 0.3, = 0.7, /i = 0.7). We use three 
difi'erent methods: a grid-based Cloud-In-Cell (GIG) 
smoothing plus a Fo urier-space estimation of the di- 
vergence (Section 3.1); an adaptive smoothing using a 



Smoothed Particle Hydrodynamics (SPH) kernel, com- 
bine d w ith a Fourier-space divergence calculation (Sec- 
tion 3.2 1; and a geometrical Delaunay tessellation es tima - 
tion ot both the density and the divergence (Section|3.3|. 



We find that each of these methods achieves differing 
degrees of success based on the nature of the length- 
scales involved: the GIG method is an Eulerian, mass- 
weighted scheme; the Delaunay tessellation is a La- 
grangian, volume-weighted scheme; and the adaptive 
mesh method is a hybrid, featuring a Lagrangian ker- 
nel with interpolation onto an Eulerian grid. We discuss 
each in turn below. 

3.1. Cloud-In-Cell and FFT 

The GIG method smooths an arbitrary distribution 
of particles onto a regularly-spaced grid by placing a 
"cloud" having the volume of a grid cell around each 
particle, so that a particle contributes to the average of 
multiple nearby cells based on the fraction of its volume 
contained in these cells. Using this weighting scheme, we 
calculate the density and displacement fields in real space 
for both 64'^ and 128'^ cells, where the displacement 
is the final {z = 0) minus initial {z = 49) particle posi- 
tions. In the case of 128^^ cells at 2; = 0, there is only one 
cell containing zero particles; we set its value to the av- 
erage value in the surrounding cells. To avoid cells with 
zero particles, this is the finest grid we use, given the 
resolution of the simulation. The divergence of the dis- 
placement can then easily be calculated in Fourier space, 
where the derivatives become simple multiples of k. 

We note here that BAG reconstruction methods also 
calculate the linear Zel'dovich displacement in Fourier 
space, which they multiply by a smoothing function S{k): 
^k.\\n = * k(5fe S{k)/k'^; for example, a Gaussian smooth- 
ing function would take the form S{k) = e^'^ ^ We 
do not apply an extra smoothing function in Fourier 
space, but note that the size of the grid cells gives the ef- 
fective smoothing length of the both the density and dis- 
placement fields. In our simulations, a 64^ cell grid has 
a length of 3.1 h^^Mpc and a 128^^ cell grid has length 
1.6 h~^Mpc, while the reconstruction methods have R — 
5, 10, an d 20 h~^Mpc, with R=10 h^^Mpc per forming 
the best ( [Eisenstein et al.||2007[ [Noh et al.||2009| ). 

In Figure [T] we show a two-dimensional histogram of 
V • \& versus S and ln(l + 5) at four different redshifts, 
for a cell size of 1.6 h~^Mpc and with a slope of -1 plot- 
ted for reference. The color scale is logarithmic so that 
the outlying cells can be seen, though the furthest out- 
liers (especially at low-z) extend beyond the range of 
the plot. At 2: = 7, the nonlinear clustering of matter 
has only begun for the initially highest-density peaks, so 
both the linear and logarithmic relations have slopes near 
-1 and are very tight. At z = 3 there is more scatter, and 
the linear relation begins to deviate very slightly from a 
slope of -1 for the majority of cells (shown in black in the 
plot) and to develop a high S tail. This trend continues 
at 2; = 1, with increased scatter in both the logarithmic 
and linear relations, and with the linear relation devi- 
ating from a slope of -1. Finally, by z = the linear 
relation deviates drastically from a slope of -1 while the 
logarithmic relation does very well, though with much 
scatter. 
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Figure 1. The divergence of the Lagrangian displacement as a 
function of S (left) and ln(l + <5) (right), calculated using CIC 
smoothing for a cell size of 1.6 h~^Mpc at 4 different redshifts. 
A slope of -1 through the origin is plotted for reference. The color 
scale is logarithmic and gives the number of CIC cells in each his- 
togram bin. 



It is not unexpected that the linear relation does not 
hold, since 1.6 h~^Mpc is by no means in the linear 
regime, but the fact that th e logarithmic relation holds 



so well echoes the results of Neyrinck et al 



and 

others that the logarithmic transform of the density field 
increases the information content of the matter power 
spectrum in the weakly nonlinear regime. When the 
CIC cell size is increased to 12.5 h^^Mpc, both the lin- 
ear and logarithmic relations become very tight, though 
the logarithmic approximation continues to do slightly 
better. Figure [2] is the same as Figure [l] but at z = 0, 
showing the CIC method for cell sizes of 1.6 h~^Mpc 
and 12.5 h~^Mpc, along with the two other methods dis- 
cussed in the next sections. 

3.2. Adaptive SPH Smoothing and FFT 

We also try an adaptiv e smo othing technique, de- 
scribed by Colombi et al. (20071, on both the density 
and the displacement helds. We assume that each parti- 
cle p traces some physical quantity Ap. For the purpose 
of this paper, Ap will be either a particle count or the 



Figure 2. The divergence of the Lagrangian displacement as a 
function of S (left) and In (1 -|- 5) (right) at 2 = 0, calculated using 
CIC smoothing for a cell size of 12.5 h~^Mpc (top) and 1.6 h~^Mpc 
(upper middle), adaptive SPH-like smoothing with a cell size of 
1.6 h~^Mpc (lower middle), and the Delaunay Tessellation Field 
Estimator (DTFE) method (bottom). A slope of -1 through the 
origin is plotted for reference. The color scale is logarithmic and 
gives the number of cells (CIC and SPH) or particles (DTFE) in 
each histogram bin. 

displacement vector 'S' between the final and initial par- 
ticle positions. We sample the adaptively-smoothed field 
onto a three-dimensional grid, with grid sites denoted 
by the integers {i,j,k). At each grid site we define the 
smoothing radius R{i, j, k) which corresponds to half the 
distance from the grid site to the M-th particle, where 
M is defined as: 



M{i,j, k) = max{M„ 



binned 



{i.j.m (7) 



and -Mbinncd(*, J, k) is the num ber of particles binned on 



the grid site (i, j, k). Following Colombi et al. (2007), we 
conservatively set -Mmin = 32, corresponding toTEe typ- 
ical number of a particle's closest neighbors. This choice 
of M{i,j,k) ensures that we always take into account a 
sufficient number of particles for making the adaptively- 
smoothed field on the grid. 

The smoothed-interpolated value A{i,j, k) at the grid 
site (i, J, fc), which corresponds to a smooth-averaging of 
the quantities Ap carried out by individual particles p, is 
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defined as: 



N 



Rii,j,kr 



Y^ApWpS 



p=l 



R{iJ,k) 



(8) 



where Wp are additional weights to be determined for 
each particle and S is the smoothing kernel function. 
We determine the weights by enforcing that the sum of 
^ij-,fc over all grid sites must be the same as the sum of 
Ap over all particles. This yields, for N particles, 



1 



N 



i,j,k i,j,k ^ ' p=l 

N N 



Rii,j,k) 



J2Ap = J2ApWpJ2 

p=l p=l i,j,k 

The last equality is true if 



1 



R{i,j,kf \R{t,j,k) 



dr, 



E 



1 



i?(i,j-fc)3 \R{t,j,k) 



(9) 



(10) 



Finally, we define the smoothi ng kernel to be th e typical 
SPH compact-support spline ( Monaghaii||1992 ), 



1 



S{x) 



3 ™2 I 3^3 
2 ' 4' 

1(2 ~x)^ l<x<2 
X > 2 



(11) 



We now compute both the mass density fluctuation 
field S{i,j, k) and the displacement field k) in Eu- 

lerian coordinates. In the case of A^-body simulations, 
the mass density fluctuation field is obtained by putting 
a mass nip = 1 on each particle p, thus 



S{i,j,k) = 



■fn{i,j,k) 
N/N^ 



1, 



(12) 



where we used the conservation of the quantity J2 ™p = 
^ rh{i,j, k), and defined Ng the number of grid sites. In 
the case of observations, one could in principle use any 
tracer of the mass that is available, such as the luminos- 
ity or the number of galaxies. The displacement field is 
obtained by averaging the displacement carried by each 
of the tracer particles. We thus obtain 



*(i,j,/c) 



'^{i,J,k) 
Hh3,k) 



(13) 



where 1 is the result of applying the adaptive filter 
to compute the number of tracers on the grid site 
{i,j,k). We note that even though in A^-body simula- 
tions l(j, j, k) = rh(i,j, k), this is not the case for obser- 
vations. 

We obtain the divergence of the displacement field in 
Fourier space, as with the previous CIC method. Fig- 
ure [2] shows V • 'S' as a function of S and ln(l + S) using 
the adaptive SPH smoothing method with a cell size of 
1.6 h~^Mpc (lower middle plot). We note that both the 
CIC and SPH methods produce very similar results. The 
additional advantage of the SPH method is that grid cells 
with zero particles are not an issue, as they are for CIC, 



since Mmin = 32 particles are always included in the bin- 
ning. The vertical artifact near ln(l -f S) — (ln(l + S)) =2 
marks the transition where Mmin becomes important, 
above which the adaptive method behaves effectively the 
same as the CIC method. Note however that this is a 
very small feature that is exaggerated by the logarithmic 
color scaling on the plot. 

3.3. Delaunay Tessellation 

The Voronoi tessellation and its dual, the Delaunay 
tessellation, divide an iV-dimensional space into a fully 
volume-covering set of cel ls based on the underlyin g par- 



ticle distribution (see e.g. van de Weygaert| ( 1994 1 ) . The 
resulting distribution of cells is such that the volume of 
the cells is inversely proportional to the density of the 
particle field at the location of each sampled particle. 



Delaunay Tessellation Field Estimator, DTFE; |Schaap 


&: van de Weygaert||2000l 


Pelupessy et al.||2003| |van de 


Weygaert & Schaap| 2009 


). In the case of the Delau- 



tetrahedra of which it makes up one of the vertices. The 
volume of one tetrahedron, T, is a function of the po- 
sitions of the three other vertices for a particle at the 
origin: 

Ft = ^xi • (x2 X X3). (14) 

The density is then given by the inverse of the volume 
of the set of tetrahedra surrounding each particle, V — 

We choose to normalize the inverse volumes by 1/{V) 
instead of (1/V^) (which gives a mean-zero density field) 
because the DTFE density distribution has a prominent 
tail at high densities. This can be seen in Figure [31 
in which we plot the distribution functions of the CIC 
ln(l+^)-(ln(l-h5)) and the DTFE ln((F)/y) = ln(l+(5). 
The CIC distribution is much more symmetrical than 
the DTFE, which has a shoulder at high densities. Note 
that the mean of the DTFE distribution has not been 
subtracted so that its peak is much closer to zero. We 
also plot a transformed DTFE distribution function in 
which we weight the original distribution function by the 
volume: P{x) = {V/{V))P{\n{{V)/V)): this effectively 
transforms the (mass-weighted) Lagrangian DTFE dis- 
tribution into a (volume-weighted) Eulerian distribution 
that closely resembles the CIC distribution. 

The Delaunay tessellation also allows a calculation of 
the divergence of ^ at the location of every particle by 
using Gauss's theorem, V • '5' dV = J^"^ ■ n dS. Re- 
placing the integrals with discrete sums, we have 

V • * ^ Vt = ^ *T • n dSr, 



(15) 



where we take to be the average of the * values of 
the three particles making the outside face of the tetra- 
hedron, and where the surface element is given by 

dSr = ^ (xi - X2) X (xi - X3) . (16) 

Thus for each particle we have a measure of both V • 'i' 
and S. 
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Figure 3. Distribution functions of ln(l + 5) — (ln(l + S)) for the 
CIC method with 1.6 h-lMpc cell sizes, \n({V) /V) for the DTFE- 
calculated volume V, and the volume-weighted DTFE transfor- 
mation (see the text for explanation). The CIC and transformed 
DTFE distributions are shifted to have zero mean, while the 
DTFE distribution is not. The good match between the CIC 
and transformed DTFE distribution functions motivate our use 
of (5 = (V) /V -1 as the DTFE density variable. 

Figure [2] shows V • \I/ as a function of 5 and ln(l + ^) at 
z = 0, calculated using the Delaunay tessellation (bot- 
tom plot) and the other methods. Both the divergence 
of ^ and the density calculated with the Delaunay tes- 
sellation lead to a large scatter at high densities beyond 
the range of the plot, with the maximum values of V • 'i' 
reaching around ±1000. This large scatter is not surpris- 
ing because neighbors in high-density regions are unlikely 
to be true Lagrangian neighbors, so their displacements 
should be fairly random; additionally, the length-scales 
associated with the highest densities are smaller than 
the scale over which we would expect the logarithmic 
approximation to hold. While the low-density success of 
the logarithmic approximation is encouraging, it seems 
that a Lagrangian kernel method is not the most appro- 
priate choice for relating the density and displacement 
fields. 

4. DISCUSSION 

Motivated by the success of the logarithmic density 
field in improving the statistics of large scale structure 
into the nonlinear regime, we have compared the linear 
relation between the density and displacement fields to a 
logarithmic approach. We have used different methods of 
measuring the divergence of the displacement field, V- ^, 
and the density field in cosmological iV-body simulations, 
and have investigated their relation for different redshifts 
and smoothing scales. The logarithmic approach fares 
better even on scales as small as ~2 h~^Mpc, where 
the linear approximation is no longer valid. Smoothing 
over a uniform grid in Eulerian coordinates seems to out- 
perform a tessellation-based Lagrangian method, though 
both work well at low to average densities. 

Because a key advantage of the logarithmic transfor- 
mation is the ability to use the density field on weakly 
nonlinear scales, this method would be most effective for 



a low-redshift galaxy survey (such that the effect of non- 
linear evolution is strongest) with a high galaxy density 
(such that the mean separation of galaxies is small and 
in the nonlinear regime). Using a logarithmic transform 
of the density to reconstruct the displacement field is 
simpler t han advanced algori thms such as MAK recon- 
struction ( Brenier et al.|2003 ) and thus may require only 
small changes to existing analysis codes, providing an 
improvement over the linear method while maintaining 
a high degree of computational efficiency. Our full in- 
vestigation of a reconstruction algorithm that uses a log- 
arithmic transform of the density field is left to future 
work and will be applied to a large suite of cosmological 
simulations. 

When applying this method to an observed galaxy dis- 
tribution, additional considerations must be accounted 
for, such as discreteness noise, galaxy bias, and redshift- 
space distortions. A log-transform can enhance the ef- 
fective shot noise caused by particle (or galaxy) discrete- 
ness, but the log-density can still be useful ly estimated 
in t he pre sence of (low-enough) shot noise. Neyrinck et 
al. (2011) found that even with a rather sparse galaxy 
sample, added power-spectrum Fisher information in the 
non-linear regime can be tapped with a log (or similar) 
transform, suggesting that shot noise can be overcome in 
the present case, as well. However, it is likely that higher 
sampling for the log-density field than for the overden- 
sity field will be necessary to resolve the fields at a given 
resolution. 

Though a full reconstruction algorithm using a log- 
density field remains to be tested, we doubt that a loga- 
rithmic transform of the density would exacerbate prob- 
lems caused by galaxy bias or redshift-space distortions, 
and these problems seem to be manageable in the case 
of BAO reconstruction. For example, it has been shown 
that linear reconstruction of the BAO peak reduces the 
effect of bias, including the induced shift in the loca- 
tion oO|il^AQ_2eak^^^lF^[201^ Addition- 



ally, Eisenstein et al. 



(|2007 ) showed that reconstruction 

of the BAO peak is eitective in both real- and redshift- 
space. However, we anticipate that Fingers of God will at 
some point prevent the push to smaller scales in redshift- 
space using a logarithmic transform. 

We have shown that a simple logarithmic transforma- 
tion of the density field provides a substantial improve- 
ment of the linear theory relation between the density 
and displacement fields on sub-linear scales. This im- 
provement is greatest at low redshift and for Eulerian 
smoothing on small (~2 h~^Mpc) scales. This result is 
in line with the finding that the logarithmic density field 
restores some lost in formation to the pow er spectrum of 
density fluctuations ( [Neyrinck et al.|[2009 |, and provides 
more evidence that a logarithmic transforin of the density 
field is a powerful density variable on weakly nonlinear 
scales. 
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